1. R requirements

1.1 Install and load packages

Load the following packages needed for analysis:

packages <- c("reshape2", "stringr", "Seurat", "tibble", "ggplot2", "dplyr", "harmony", "clusterProfiler", "ComplexHeatmap", "ggpubr", "ggrepel", "ggnewscale", "Matrix", "useful", "RColorBrewer", "org.Hs.eg.db", "tidyr", "circlize", "patchwork", "pheatmap", "gridExtra", "grid", "cowplot", "viridis", "fmsb")

lapply(packages, require, character.only = TRUE)

1.2 Load functions

Functions for downstream analysis

source("Functions.R")
hallmark_genes_symbols <- read.gmt("./files/h.all.v7.4.symbols.gmt")
hallmark_genes         <- read.gmt("./files/h.all.v7.4.entrez.gmt")

reactome_genes_symbols <- read.gmt("./files/c2.cp.reactome.v7.4.symbols.gmt")
reactome_genes         <- read.gmt("./files/c2.cp.reactome.v7.4.entrez.gmt")

KEGG_genes_symbols<-read.gmt("./files/c2.cp.kegg.v7.4.symbols.gmt")
KEGG_genes<-read.gmt("./files/c2.cp.kegg.v7.4.entrez.gmt")

gmtfile_hallmarks <- hallmark_genes

1.3 Load colors

Used colors in this manuscript.

color_clusters<-c(brewer.pal(n = 9,name = "Set1"),
                  brewer.pal(n = 8,name = "Set2"),
                  brewer.pal(n = 12,name = "Set3"),
                  brewer.pal(n = 12,name = "Paired"),
                  brewer.pal(n = 9,name = "Pastel1"),
                  brewer.pal(n = 8,name = "Pastel2"),
                  brewer.pal(n = 8,name = "Accent"))

colors_celltype <- c("CD8 T cells" = "#B02457",
                  "CD4 T cells" = '#ea7800',
                  "NK cells" = '#66357D',
                  "Megakaryocytes" = '#5a6b22',
                  "B cells" = "#fcc113",
                  "Monocytes" = "#9dcfcf",
                  "Plasmablasts" = '#17315f',
                  "mDCs" = "#C19BA8",
                  "pDCs" = "#BF40A8",
                  "Erythrocytes" = "red",
                  "Prol. cells" = '#9eceaa')

colors_monocytes <- c(  "IFIhi" = "#6d3b8f",  
                        "IL1Bhi" = "#b50e3b", 
                        "S100Ahi" = "#4D39D1",
                        "CXCLhi" = "#fcc113",  
                        "Dexa_response" = "#41b6c4",
                        "classical" = "#FBB99C",
                        "THBDhi" = "#e04b07",
                        "PF4+" = "#e54585",
                        "HBhi" = "lightgrey",
                        "non_classical" = "#FF99FF")

colors_severity <- c("severe" = "#CD5C5C",
                     "moderate" = "#5B9BD5")

colors_severity_ctrl <- c("severe_ctrl" = "indianred",
                     "moderate_ctrl" = "orange")

2. Load Datasets

2.1 PBMC data

seurat_PBMC <- readRDS("COVID_Dexa_seurat_PBMC.RDS")
seurat_PBMC
## An object of class Seurat 
## 25498 features across 114181 samples within 2 assays 
## Active assay: RNA (25486 features, 2000 variable features)
##  1 other assay present: HTO
##  2 dimensional reductions calculated: pca, umap

2.2 monocyte data

seurat_monocytes <- readRDS("COVID_Dexa_seurat_monocytes.RDS")
seurat_monocytes
## An object of class Seurat 
## 25498 features across 23648 samples within 2 assays 
## Active assay: RNA (25486 features, 1000 variable features)
##  1 other assay present: HTO
##  3 dimensional reductions calculated: pca, umap, harmony

3. Figure 1

3.1 UMAP by celltype

Visualization of cell types on the UMAP embedding.

seurat_PBMC$cluster_labels_res.0.4 <- factor(seurat_PBMC$cluster_labels_res.0.4, levels = c( "Monocytes", "mDCs", "pDCs", 
                                                                                             "Megakaryocytes", "CD4 T cells", 
                                                                                             "CD8 T cells", "NK cells", "B cells",
                                                                                             "Plasmablasts", "Prol. cells",  "Erythrocytes"))

DimPlot(seurat_PBMC, reduction = "umap", group.by = "cluster_labels_res.0.4", cols = colors_celltype, raster = F)

3.2 UMAP by severity

Subsetting the data to controls only (untreated COVID-19) and highlighting severity.

seurat_ctrls <- subset(seurat_PBMC, group %in% c("moderate_ctrl", "severe_ctrl"))

p1 <-DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(14000), sizes.highlight = 0.1, cols.highlight = "#FFB651", cols = "lightgrey")+NoLegend()+ggtitle("moderate")
p2 <- DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(row.names(seurat_ctrls@meta.data[seurat_ctrls@meta.data$group == "severe_ctrl",]), 14000), sizes.highlight = 0.1, cols.highlight = "#CD5C5C", cols = "lightgrey")+NoLegend()+ggtitle("severe")

p1 + p2

3.3 UMAP by early/late

Subsetting the data to controls only (untreated COVID-19) and highlighting timepoints.

seurat_ctrls$timepoints <- ifelse(seurat_ctrls$SO.to.sample > 10, "late", "early" )

p1 <- DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(row.names(seurat_ctrls@meta.data[seurat_ctrls@meta.data$timepoints == "early",]), 17000), sizes.highlight = 0.1, cols.highlight = "#c51b7d", cols = "lightgrey", order = F, raster = F)+NoLegend()+ggtitle("Early  <= 10d")
p2 <-DimPlot(seurat_ctrls, reduction = "umap", cells.highlight = sample(row.names(seurat_ctrls@meta.data[seurat_ctrls@meta.data$timepoints == "late",]), 17000), sizes.highlight = 0.1, cols.highlight = "#92c5de", cols = "lightgrey", order = F, raster = F)+NoLegend()+ggtitle("Late >10d")

p1+p2

3.4 Marker by cell type

Top marker by cell types.

seurat_PBMC$cluster_labels_res.0.4 <- factor(seurat_PBMC$cluster_labels_res.0.4, levels = c( "Monocytes", "mDCs", "pDCs", 
                                                                                             "Megakaryocytes", "CD4 T cells", 
                                                                                             "CD8 T cells", "NK cells", "B cells",
                                                                                             "Plasmablasts", "Prol. cells",  "Erythrocytes"))

Idents(seurat_PBMC) <- "cluster_labels_res.0.4"

cluster.markers.wilcox_cluster_labels_res.0.4 <- FindAllMarkers(object = seurat_PBMC, 
                                         only.pos = TRUE, 
                                         min.pct = 0.2, 
                                         logfc.threshold = 0.5,
                                         min.diff.pct = 0.2,
                                         test.use = "wilcox", max.cells.per.ident = 7000, verbose = F
                                         )

top <- cluster.markers.wilcox_cluster_labels_res.0.4 %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

DotPlot(seurat_PBMC, features = rev(unique(top$gene)), 
              group.by = "cluster_labels_res.0.4", col.max = 1.5, col.min = -1.5)+ 
              scale_color_gradientn(colors = rev(RColorBrewer::brewer.pal(n =20, name = "RdBu")))+
              theme(axis.text.x = element_text(size = 14, angle = 90, hjust = 1, vjust =0.5), 
                    legend.position = "bottom", axis.text.y = element_text(face = "italic", size = 12))+
              coord_flip()

3.5 Monocyte HLA and alarmin

Checking “ground truth” for HLA and alarmin expression in untreated COVID-19 patients monocytes.

seurat_ctrls_monos <- subset(seurat_monocytes, group %in% c("severe_ctrl", "moderate_ctrl"))

# selected HLA 
p1 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("HLA-DRA"), pt.size = 0, cols = colors_severity_ctrl)+  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+theme( axis.title.y = element_blank(), axis.title.x = element_blank())
p2 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("HLA-DRB1"), pt.size = 0, cols = colors_severity_ctrl)+  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme( axis.title.y = element_blank(), axis.title.x = element_blank())
p3 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("HLA-DPA1"), pt.size = 0, cols = colors_severity_ctrl)+  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())

# selected alarmins
p4 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("S100A8"), pt.size = 0, cols = colors_severity_ctrl)+  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p5 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("S100A9"), pt.size = 0, cols = colors_severity_ctrl)+  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p6 <- VlnPlot(seurat_ctrls_monos, group.by = "group", features = c("S100A12"), pt.size = 0, cols = colors_severity_ctrl)+  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+NoLegend()+NoLegend()+theme(axis.title.y = element_blank(), axis.title.x = element_blank())


p1+p2+p3+p4+p5+p6+plot_layout(ncol = 6)

Testing significance for selected genes

df <- FindMarkers(seurat_ctrls_monos, group.by = "group", ident.1 = "severe_ctrl", ident.2 = "moderate_ctrl",  features = c("HLA-DRA", "HLA-DRB1", "HLA-DPA1", "S100A8", "S100A9", "S100A12"))
df$gene <- row.names(df)

df[sort(df$gene), c("avg_log2FC", "p_val_adj")]
##          avg_log2FC     p_val_adj
## HLA-DPA1 -0.9061705 6.252979e-208
## HLA-DRA  -0.9458845 2.937397e-278
## HLA-DRB1 -1.3566956  0.000000e+00
## S100A12   1.1022834 3.743919e-301
## S100A8    1.0036270  0.000000e+00
## S100A9    0.5324372 2.161001e-159

3.6 IFN module monocytes

IFN-gene expression enrichment in monocytes by timepoints.

seurat <- subset(seurat_monocytes, group_severity %in% c("moderate", "severe"))
seurat$order_group <- paste(seurat$order, seurat$group, sep = "_")

seurat$order_group <- gsub(seurat$order_group, pattern = "early_severe_Dexa", replacement = "early")
seurat$order_group <- gsub(seurat$order_group, pattern = "early_severe_ctrl", replacement = "early")
seurat$order_group <- gsub(seurat$order_group, pattern = "early_moderate_Dexa", replacement = "early")
seurat$order_group <- gsub(seurat$order_group, pattern = "early_moderate_ctrl", replacement = "early")

seurat$order_group <- gsub(seurat$order_group, pattern = "late_severe_Dexa", replacement = "late_Dexa")
seurat$order_group <- gsub(seurat$order_group, pattern = "late_moderate_Dexa", replacement = "late_Dexa")

seurat$order_group <- gsub(seurat$order_group, pattern = "late_severe_ctrl", replacement = "late_Ctrl")
seurat$order_group <- gsub(seurat$order_group, pattern = "late_moderate_ctrl", replacement = "late_Ctrl")

seurat$sev_order_group <- paste(seurat$group_severity, seurat$order_group, sep = "_")

IFNA_signature <- c("IFI6", "ISG15", "IFITM1", "ISG20", 
                    "IFI27", "IFI30", "IFIH1", "IFIT1", 
                    "IFIT2", "IFIT3", "IFITM2", "IFITM3", 
                    "XAF1", "MX1", "MX2")

seurat <- AddModuleScore(seurat, features = list(IFNA_signature), name = "IFNA_signature")

df <- seurat@meta.data[,c("donorID", "sampleID", "sev_order_group", "IFNA_signature1")]
df %>% group_by(sev_order_group) %>% summarize(mean = mean(IFNA_signature1)) -> df.mean
df <- merge(df.mean, df, by = "sev_order_group")

p <- ggplot(df, aes(x = sev_order_group, y = IFNA_signature1))+
  geom_violin(aes(fill = mean))+
  stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0)+
  theme_classic()+
  theme(panel.grid.major = element_line("white"),
        panel.grid.minor = element_line("white"),
        axis.title.y = element_text( size = 12),
        axis.text.y = element_text( size = 12),
        axis.title.x = element_text( size = 12),
        axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
        legend.text =  element_text( size = 12),
        plot.title = element_text(size = 14, hjust = 0.5))+
  xlab("")+
  ylab("Module Score")+
  ggtitle(label = "IFNA_signature", subtitle = paste("n genes = ", length(IFNA_signature)))+
  NoLegend()

p

3.7 PBMC DE call

DEG call for Dexa vs. control for both moderate and severe COVID-19.

seurat <- subset(seurat_PBMC, order == "late" & group_severity %in% c("moderate", "severe"))

# define comparisons
comparisons<-data.frame(comparison=c("sev_Dex_vs_ctrl","mod_Dex_vs_ctrl"),
               group1=c("severe_Dexa","moderate_Dexa"),
               group2=c("severe_ctrl","moderate_ctrl"))

#define cell types to include
meta<-seurat@meta.data
celltype_to_show <- c("Monocytes", "CD4 T cells", "CD8 T cells", "NK cells",  "B cells")

DEgenes.all<-data.frame()

for(i in comparisons$comparison){
  groups_oi<-as.vector(comparisons[comparisons$comparison==i,])
  seurat_tmp<-subset(seurat, subset = group %in% groups_oi)

  for(j in celltype_to_show){
    Idents(seurat_tmp)<-"cluster_labels_res.0.4"
    seurat_celltype<-subset(seurat_tmp, subset = cluster_labels_res.0.4==j)
    Idents(seurat_celltype)<-"group"
    tmp <-FindMarkers(object = seurat_celltype,
                      ident.1 = comparisons[comparisons$comparison==i,]$group1,
                      ident.2 = comparisons[comparisons$comparison==i,]$group2,
                      only.pos = FALSE,
                      min.pct = 0.1,
                      logfc.threshold = 0.25,
                      test.use = "wilcox")
    #add information
    tmp$celltype<-j
    tmp$cluster<-paste0(j,"_",i)
    tmp$comparison<-i
    tmp$celltype_comparison<-paste0(j,"_",i)
    tmp$gene<-row.names(tmp)
    tmp$regulation<-ifelse(tmp$avg_log2FC>0,"up", "down")
    tmp$significance<-ifelse(tmp$p_val_adj<0.05,"<0.05", "ns")
 
    DEgenes.all<-rbind(DEgenes.all,tmp)
  }
}

rm(seurat_tmp, seurat_celltype,groups_oi, tmp)
DEgenes<-DEgenes.all[DEgenes.all$p_val_adj<0.05,]

DEgenes$celltype <- factor(DEgenes$celltype, levels = c("Monocytes","B cells", "CD4 T cells", "CD8 T cells",  "NK cells"))

p <- ggplot(DEgenes, aes(comparison, fill=regulation))+
  geom_bar(position = position_dodge(), color = "black")+
  scale_fill_manual(values = c(up = "#F59F9F", down = "#76B2ED"))+
  scale_alpha_manual(values = c("<0.05"=1, "ns"=0.8))+
  theme_linedraw()+
  geom_text(stat='count',
              aes(label=..count..),
              position=position_dodge2(width=0.8), size=3, vjust = -0.3, hjust = 0.6)+
  theme(axis.text.x = element_text(angle = 90, hjust=1,vjust = 0))+
  ggtitle("DE genes (p.adj<0.05)")+
  xlab("")+
  ylim(0,260)+
  facet_wrap(.~celltype,nrow = 1)+
  theme(strip.background =element_rect(fill="grey"))+
  theme(strip.text = element_text(color = "black",size = 10))

p

3.8 Common DEG severe

Identification of genes that are DEG in at least 3 cell types in severe COVID-19 Dexa vs. control.

DEgenes_sev <- DEgenes[DEgenes$comparison == "sev_Dex_vs_ctrl", ]

# check commond DEG down
DEgenes_down <- DEgenes_sev[DEgenes_sev$regulation == "down", ]

tmp_all_down <- c()
for(i in unique(DEgenes_down$celltype)){
  tmp <- DEgenes_down[DEgenes_down$celltype == i, ]
  name <- gsub(pattern = "_sev_Dex_vs_ctrl_down", replacement = "", x = tmp$cluster_regulation)
  name <- unique(name)
  name <- gsub(name, pattern = " ", replacement = "_")
  tmp <- tmp$gene
  tmp_all_down <- c(tmp_all_down, tmp)
}
df <- as.data.frame(tmp_all_down)
colnames(df) <- "genes"

genes_multiple_celltypes_down <- vector()

for(i in unique(df$genes)){
  tmp <- df[df$genes == i, ]
  tmp <- as.data.frame(tmp)
  if(nrow(tmp) >= 3){
    genes_multiple_celltypes_down <- c(i, genes_multiple_celltypes_down)
  } 
}

# check common DEG down
DEgenes_sev <- DEgenes[DEgenes$comparison == "sev_Dex_vs_ctrl", ]
DEgenes_up <- DEgenes_sev[DEgenes_sev$regulation == "up", ]

tmp_all_up <- c()
for(i in unique(DEgenes_up$celltype)){
  tmp <- DEgenes_up[DEgenes_up$celltype == i, ]
  name <- gsub(pattern = "_sev_Dex_vs_ctrl_up", replacement = "", x = tmp$cluster_regulation)
  name <- unique(name)
  name <- gsub(name, pattern = " ", replacement = "_")
  tmp <- tmp$gene
  tmp_all_up <- c(tmp_all_up, tmp)
}
df <- as.data.frame(tmp_all_up)
colnames(df) <- "genes"

genes_multiple_celltypes_up <- vector()
for(i in unique(df$genes)){
  tmp <- df[df$genes == i, ]
  tmp <- as.data.frame(tmp)
  if(nrow(tmp) >= 3){
    genes_multiple_celltypes_up <- c(i, genes_multiple_celltypes_up)
  } 
}


list(up = genes_multiple_celltypes_up,
     down = genes_multiple_celltypes_down)
## $up
## [1] "HLA-F"   "SF3A2"   "RPS5"    "RPL39"   "TXNIP"   "TSC22D3"
## 
## $down
##  [1] "ZNF331"     "NR4A2"      "EEF1A1"     "SNORD3A"    "EIF5A"     
##  [6] "AL590867.2" "HMGB2"      "RPL9P9"     "IFITM1"     "MT-ATP8"   
## [11] "RPL10"      "RPL6"       "PTMA"       "RPL4"       "RPL14"     
## [16] "MTND2P28"   "LITAF"      "AP005202.1" "AC109326.1" "AL365357.1"
## [21] "FTH1"

common DEG severe include TSC22D3 and TXNIP (up) and FTH1 and IFITM1 (down).

seurat <- subset(seurat_PBMC, order == "late" & group_severity %in% c("severe"))

seurat <- subset(seurat, cluster_labels_res.0.4 %in% c("Monocytes","B cells", "CD4 T cells", 
                                                       "CD8 T cells",  "NK cells"))

seurat$cluster_labels_res.0.4 <- factor(seurat$cluster_labels_res.0.4 , levels = c("Monocytes","B cells", "CD4 T cells", 
                                                                                   "CD8 T cells",  "NK cells"))

p1 <- VlnPlot(seurat, features = "TSC22D3", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
  stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p2 <- VlnPlot(seurat, features = "TXNIP", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
  stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p3 <- VlnPlot(seurat, features = "FTH1", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
  stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
p4 <- VlnPlot(seurat, features = "IFITM1", group.by = "cluster_labels_res.0.4", split.by = "treatment", split.plot = T, pt.size = 0, cols = c("#984EA3", "#5AAE61"))+
  stat_summary(fun.y = mean, geom='point', size = 3, colour = "black", shape = 95, position = position_dodge(0.75))+NoLegend()+theme(axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

p1+p2+p3+p4+plot_layout(ncol = 4)

statistics check- to identidy DEG.

genes <- c("TSC22D3", "TXNIP", "FTH1", "IFITM1")
tmp <- DEgenes_sev[DEgenes_sev$gene %in% genes, ]
tmp <- tmp[order(tmp$gene, decreasing = F), ]
tmp[, c("celltype", "gene", "p_val_adj")]
##             celltype    gene     p_val_adj
## FTH1       Monocytes    FTH1 6.317354e-109
## FTH11    CD4 T cells    FTH1  5.448513e-34
## FTH12    CD8 T cells    FTH1  4.043850e-07
## FTH13       NK cells    FTH1  2.339069e-13
## FTH14        B cells    FTH1  4.965424e-37
## IFITM1   CD4 T cells  IFITM1  8.021860e-35
## IFITM11  CD8 T cells  IFITM1  1.947397e-07
## IFITM12     NK cells  IFITM1  2.496873e-10
## TSC22D3    Monocytes TSC22D3 5.211175e-130
## TSC22D31 CD4 T cells TSC22D3  5.756550e-55
## TSC22D32 CD8 T cells TSC22D3  2.721599e-06
## TSC22D33    NK cells TSC22D3  6.003117e-03
## TSC22D34     B cells TSC22D3  4.622287e-39
## TXNIP      Monocytes   TXNIP  2.363646e-41
## TXNIP1   CD8 T cells   TXNIP  1.668797e-09
## TXNIP2      NK cells   TXNIP  3.626118e-05
## TXNIP3       B cells   TXNIP  2.128502e-42

3.9 GOEA PBMC

3.9.1 Common GOEA in multiple cell types

DEgenes$cluster_regulation<-paste(DEgenes$celltype_comparison, DEgenes$regulation, sep="_")
DE.list<-list()
for(i in unique(DEgenes$cluster_regulation)){
  DE.list[[i]]<-DEgenes[DEgenes$cluster_regulation==i,]$gene
}

GOEA_celltypes_oi<-GOEA.of.list(seurat_object= seurat, 
                DE.list,
             GeneSets =c("GO","Hallmark", "KEGG",  "Reactome"
                         ),
             GOntology = "BP",
             pCorrection = "bonferroni",      # choose the p-value adjustment method (Steffi: bonferroni)
             pvalueCutoff = 0.05,      # set the unadj. or adj. p-value cutoff (depending on correction method)
             qvalueCutoff = 0.05       # set the q-value cutoff (FDR corrected)
)

GOEA_celltypes_oi$KEGGresults <- NULL


GOEA_celltypes_oi$HALLMARKresults$database <- "Hallmark"
GOEA_celltypes_oi$GOresults$database <- "GO"
GOEA_celltypes_oi$REACTOMEresults$database <- "Reactome"

GOEA_celltypes_oi$REACTOMEresults$Description <- gsub(GOEA_celltypes_oi$REACTOMEresults$Description, pattern = "REACTOME_", replacement = "")

GOEA_all <- rbind(GOEA_celltypes_oi$GOresults, 
                  rbind(GOEA_celltypes_oi$REACTOMEresults, GOEA_celltypes_oi$HALLMARKresults))

GOEA_all$geneSet <- row.names(GOEA_all)

GOEA_all$group <-  sapply(strsplit(GOEA_all$geneSet, split = "_"), `[`, 2)


GOEA_all %>% group_by(Condition) %>% arrange(p.adjust, .by_group = TRUE) -> GOEA_all
GOEA_all %>% group_by(Condition) %>% arrange(Count, .by_group = TRUE) %>% top_n(n= 10, wt = p.adjust) %>% pull(Description) ->terms
GOEA_all <- GOEA_all[GOEA_all$Description %in% terms,]
GOEA_all$Description <- ifelse(nchar(GOEA_all$Description) > 50,
                        paste(substr(GOEA_all$Description, 1, 50), "[...]", sep=""),
                        GOEA_all$Description)
GOEA_all$Description <- factor(GOEA_all$Description, levels = rev(unique(GOEA_all$Description)))

#change color code of dots
GOEA_all$DE.regulation<-"no"
if(nrow(GOEA_all[grep("_up", GOEA_all$Condition),])>0){
  GOEA_all[grep("_up", GOEA_all$Condition),]$DE.regulation<-"up"
}
if(nrow(GOEA_all[grep("_down", GOEA_all$Condition),])>0){
  GOEA_all[grep("_down", GOEA_all$Condition),]$DE.regulation<-"down"
}


GOEA_all$Cond_short <- gsub(GOEA_all$Condition, pattern = "_.*", replacement = "")



GOEA_all_small <- GOEA_all[, c("Cond_short", "Description", "Count", "group", "DE.regulation")]
nrow(GOEA_all_small)
## [1] 302
df <- GOEA_all_small %>% group_by(Description, group, DE.regulation) %>% summarize(count = n())

GOEA_all_small <- GOEA_all_small[GOEA_all_small$Description %in% df[df$count >= 3, ]$Description,  ]

GOEA_all_small <- GOEA_all_small[GOEA_all_small$Cond_short %in% c("Monocytes", "B cells", "CD4 T cells", "CD8 T cells", "NK cells"), ]

GOEA_all_small$Cond_short <- factor(GOEA_all_small$Cond_short, levels = c("Monocytes", "B cells", "CD4 T cells", "CD8 T cells", "NK cells"))

p <- ggplot(data = GOEA_all_small ,aes(x=Cond_short, y=Description, size=Count)) +
  geom_point(fill = "black",pch=21)+theme_classic()+
  theme_linedraw()+
  theme(axis.text.y = element_text(color = "black"),
        axis.text.x = element_text(angle=90, vjust=0.5,hjust=1, color = "black"),
        text = element_text(size = 12))+
  scale_fill_manual(values= c(sev = "#CD5C5C", mod = "#FFB651"))+
  xlab("")+ylab("")+
  scale_radius(breaks = waiver())+
  facet_wrap(DE.regulation ~ group, scales = "fixed")

p

3.9.2 NFKB down spider

selection of HALLMARK_TNFA_SIGNALING_VIA_NFKB as the best term for down regulated genes and checking enrichment again in a spyder plot.

df <- GOEA_all_small[GOEA_all_small$Description %in% "HALLMARK_TNFA_SIGNALING_VIA_NFKB" & GOEA_all_small$DE.regulation == "down", ]

extension <- data.frame("Cond_short" = c("NK cells", "CD4 T cells", "CD8 T cells",
                                         "CD4 T cells", "NK cells"), 
                        "Description" = c("HALLMARK_TNFA_SIGNALING_VIA_NFKB",
                                          "HALLMARK_TNFA_SIGNALING_VIA_NFKB", 
                                          "HALLMARK_TNFA_SIGNALING_VIA_NFKB",
                                          "HALLMARK_TNFA_SIGNALING_VIA_NFKB", "HALLMARK_TNFA_SIGNALING_VIA_NFKB"),
                        "Count" = c(0, 0, 0, 0, 0),
                        "group" = c("mod", "mod", "mod", "sev", "sev"),
                        "DE.regulation" = c("down", "down", "down", "down", "down"))

df <- rbind(df, extension)
df$Cond_short <- factor(df$Cond_short, levels = c("Monocytes", "B cells", "CD4 T cells", "CD8 T cells", "NK cells"))
df_mod <- df[df$group == "mod", ]
df_sev <- df[df$group == "sev", ]

max <- ceiling(max(c(max(df_sev$Count), max(df_mod$Count))))
data <- matrix(rep(0, 5), ncol = 5)
colnames(data) <- c("Monocytes", "B cells", "NK cells", "CD8 T cells", "CD4 T cells")
tmp <- rbind(rep(40, 5), data,  df_mod$Count[match(colnames(data), df_mod$Cond_short)])
tmp <- as.data.frame(tmp)
tmp[is.na(tmp)] <- 0

radarchart(tmp,  axistype=1 , 
    #custom polygon
    pcol=rgb(1, 0.714, 0.318,0.9) , pfcol=rgb(1, 0.714, 0.318,0.5) , plwd=4 , 
    #custom the grid
    cglcol="darkgrey", cglty=1, axislabcol="black", cglwd=1, caxislabels=seq(0,40,10), seg = 4, title = "moderate NF-kB TF prediciton",
    #custom labels
    vlcex=1 )

data <- matrix(rep(0, 5), ncol = 5)
colnames(data) <- c("Monocytes", "B cells", "NK cells", "CD8 T cells", "CD4 T cells")
tmp <- rbind(rep(40, 5), data,  df_sev$Count[match(colnames(data), df_sev$Cond_short)])
tmp <- as.data.frame(tmp)
tmp[is.na(tmp)] <- 0
radarchart(tmp,  axistype=1 , 
    #custom polygon
    pcol=rgb(0.804, 0.361, 0.361, 0.9) , pfcol=rgb(0.804, 0.361, 0.361,0.5) , plwd=4 , 
    #custom the grid
    cglcol="darkgrey", cglty=1, axislabcol="black", cglwd=1, caxislabels=seq(0,40,10), seg = 4, title = "severate NF-kB TF prediciton",
    #custom labels
    vlcex=1 )

4. Figure 2

4.1 FC-FC plot Dexa core

To identify the Core Dexa signature in monocytes we compare DEG from the comparison Dexa vs control for moderate and severe COVID-19.

DEG_list <- list()

for(i in c("moderate", "severe")){
  
  tmp <- subset(seurat_monocytes, subset = group_severity == i)

  for (x in c("late")) {
    
    tmp_sample_order <- subset(tmp, subset = order == x)
    tmp_sample_order <- SetIdent(tmp_sample_order, value = "group")
    
      ctrl <- unique(tmp_sample_order$group)[grep(unique(tmp_sample_order$group), pattern = "ctrl")]
      
      idx <- ifelse(grep(unique(tmp_sample_order$group), pattern = "ctrl") == 1, 2, 1)
      treat <- unique(tmp_sample_order$group)[idx]
      
      tryCatch({
        DEG<-FindMarkers(tmp_sample_order,
                         ident.1 =  treat,
                         ident.2 = ctrl,
                         min.pct = 0.1,
                         logfc.threshold = 0,
                         test.use = "wilcox",
                         slot = "data",
                         only.pos = F,
                         verbose = F)
        DEG$Gene <- rownames(DEG)
        DEG$group_severity <- paste(i)
        DEG$sample_comparison <- paste(x)
        DEG$comparison <- paste0(i, "_", x, "_treatment_vs_ctrl")
        DEG_list[[paste0(i, "_", x, "_treatment_vs_ctrl")]]<-DEG
      }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
 }
}


tmp1 <- DEG_list$moderate_late_treatment_vs_ctrl[intersect(row.names(DEG_list$moderate_late_treatment_vs_ctrl), row.names(DEG_list$severe_late_treatment_vs_ctrl)),]
tmp2 <- DEG_list$severe_late_treatment_vs_ctrl[intersect(row.names(DEG_list$moderate_late_treatment_vs_ctrl), row.names(DEG_list$severe_late_treatment_vs_ctrl)),]
colnames(tmp1) <- paste(colnames(tmp1), "moderate", sep = "_")
colnames(tmp2) <- paste(colnames(tmp2), "severe", sep = "_")
data <- cbind(tmp1, tmp2)

data$Gene <- row.names(data)
data_rank_dexa_mod_severe <- data
data_rank_dexa_mod_severe$regulation <- ifelse(data_rank_dexa_mod_severe$p_val_adj_moderate < 0.05 & data_rank_dexa_mod_severe$p_val_adj_severe < 0.05, "DE_both", "DE_one")
data_rank_dexa_mod_severe$keep <- ifelse(data_rank_dexa_mod_severe$p_val_adj_moderate >= 0.05 & data_rank_dexa_mod_severe$p_val_adj_severe >= 0.05, "remove", "keep")
data_rank_dexa_mod_severe[data_rank_dexa_mod_severe$keep %in% "remove", "regulation"] <- "n.s."

df <- data_rank_dexa_mod_severe

down_label <- df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_both",]
down_label_mod <- down_label[order(down_label$avg_log2FC_moderate), "Gene"][1:15]
down_label_sev <- down_label[order(down_label$avg_log2FC_severe), "Gene"][1:15]
down_label_gene <- unique(down_label_sev, down_label_mod)

up_label <- df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_both",]
up_label_mod <- up_label[order(up_label$avg_log2FC_moderate, decreasing = T), "Gene"][1:15]
up_label_sev <- up_label[order(up_label$avg_log2FC_severe, decreasing = T), "Gene"][1:15]
up_label_gene <- unique(up_label_sev, up_label_mod)


p <- ggplot(df, aes(x = avg_log2FC_moderate, y = avg_log2FC_severe))+
            geom_point(data =df,
                       pch = 21,size = 2, color = "black", fill = "white")+
  
            geom_point(data= df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_one",],    
                       aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#c6dbef")+
            geom_point(data= df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_one",],    
                       aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#fcc5c0")+
           
  
            geom_text_repel(data= df[df$Gene %in% down_label_gene, ],
                                    aes(x=avg_log2FC_moderate, y=avg_log2FC_severe, label = Gene),
                          size = 4, color = "black", fontface = 'italic', max.overlaps = 40)+
            geom_point(data= df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_both",],    
                       aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#2171b5")+

            geom_text_repel(data= df[df$Gene %in% up_label_gene, ],
                                     aes(x=avg_log2FC_moderate, y=avg_log2FC_severe, label = Gene),
                            size = 4, color = "black", fontface = 'italic', max.overlaps = 40)+
            geom_point(data= df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_both",],                          aes(x=avg_log2FC_moderate, y=avg_log2FC_severe),size = 3, pch = 21, fill = "#cb181d")+
            theme_bw()+ 
            theme_linedraw()+
            theme(panel.grid.major = element_line("white"),
                  panel.grid.minor = element_line("white"),
                  axis.title.y = element_text( size = 12),
                  axis.text.y = element_text( size = 10),
                  axis.title.x = element_text( size = 12),
                  axis.text.x = element_text( size = 10),
                  strip.background = element_rect(fill=c("lightgrey")),
                  strip.text = element_text(colour = "black", face = "bold", size = 12), 
                  plot.title = element_text(face = "bold", size = 12))+
                  xlab("log2FC moderate")+
                  ylab("log2FC severe")+
            ggtitle("monocytes dexa: severe vs moderate")+ 
            geom_hline(yintercept = 0)+ 
            geom_vline(xintercept = 0)+
            geom_abline(linetype = 2)+
            xlim(c(-1.6,3.7))+
            ylim(c(-1.7, 2.1))

p

4.2 Dexa-induced gene expression violins

seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat_monocytes_late, group_severity %in% c("moderate", "severe"))

seurat$group <- factor(seurat$group, levels = c("moderate_ctrl", "moderate_Dexa",
                                                     "severe_ctrl", "severe_Dexa"))
p1<-VlnPlot(seurat, features = c("TSC22D3"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)
p2<-VlnPlot(seurat, features = c("IL1R2"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p3<-VlnPlot(seurat, features = c("CD163"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p4<-VlnPlot(seurat, features = c("SAP30"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p5<- VlnPlot(seurat, features = c("PER1"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())
p6<- VlnPlot(seurat, features = c("ZFP36L2"), group.by = "group", pt.size = 0, cols = c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61"))+NoLegend()+
stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+theme(axis.title.y = element_blank())

p1+p2+p3+p4+p5+p6+plot_layout(ncol = 3)

checking significance for DEG:

tmp <- DEG_list$moderate_late_treatment_vs_ctrl
tmp <- tmp[tmp$Gene %in% c("TSC22D3", "IL1R2", "CD163", "SAP30", "PER1", "ZFP36L2"), c("group_severity", "Gene", "p_val_adj", "avg_log2FC")]
tmp[order(tmp$Gene), ]
##         group_severity    Gene     p_val_adj avg_log2FC
## CD163         moderate   CD163 5.191492e-122  1.1818631
## IL1R2         moderate   IL1R2  0.000000e+00  2.6483134
## PER1          moderate    PER1  1.735679e-78  0.8029548
## SAP30         moderate   SAP30 4.544743e-135  1.2053391
## TSC22D3       moderate TSC22D3  3.557426e-91  0.7122654
## ZFP36L2       moderate ZFP36L2  7.352414e-68  0.9031032
tmp <- DEG_list$severe_late_treatment_vs_ctrl
tmp <- tmp[tmp$Gene %in% c("TSC22D3", "IL1R2", "CD163", "SAP30", "PER1", "ZFP36L2"), c("group_severity", "Gene", "p_val_adj", "avg_log2FC")]
tmp[order(tmp$Gene), ]
##         group_severity    Gene     p_val_adj avg_log2FC
## CD163           severe   CD163  2.929259e-20  0.5129026
## IL1R2           severe   IL1R2  5.249462e-87  1.7486982
## PER1            severe    PER1  1.472531e-07  0.3160355
## SAP30           severe   SAP30  5.188564e-23  0.4672367
## TSC22D3         severe TSC22D3 2.212959e-117  0.6846452
## ZFP36L2         severe ZFP36L2  2.804138e-08  0.4505909

4.3 Core Heatmap

Core Dexa DEG as heatmap (top 20 up and down).

df <- data_rank_dexa_mod_severe

down_label <- df[df$avg_log2FC_moderate <0 & df$avg_log2FC_severe <0 & df$regulation %in% "DE_both",]
down_label_mod <- down_label[order(down_label$avg_log2FC_moderate), "Gene"][1:20]
down_label_sev <- down_label[order(down_label$avg_log2FC_severe), "Gene"][1:20]
down_label_gene <- unique(down_label_sev, down_label_mod)

up_label <- df[df$avg_log2FC_moderate >0 & df$avg_log2FC_severe >0 & df$regulation %in% "DE_both",]
up_label_mod <- up_label[order(up_label$avg_log2FC_moderate, decreasing = T), "Gene"][1:20]
up_label_sev <- up_label[order(up_label$avg_log2FC_severe, decreasing = T), "Gene"][1:20]
up_label_gene <- unique(up_label_sev, up_label_mod)

df_1 <- DEG_list$moderate_late_treatment_vs_ctrl[row.names(DEG_list$moderate_late_treatment_vs_ctrl) %in% c(up_label_gene, down_label_gene),]
colnames(df_1) <- paste(colnames(df_1), "moderate", sep = "_")

df_2 <- DEG_list$severe_late_treatment_vs_ctrl[row.names(DEG_list$severe_late_treatment_vs_ctrl) %in% c(up_label_gene, down_label_gene),]
colnames(df_2) <- paste(colnames(df_2), "severe", sep = "_")

df_1 <- df_1[c(up_label_gene, down_label_gene),]
df_2 <- df_2[c(up_label_gene, down_label_gene),]

df_all <- cbind(df_1, df_2)

df_all$Gene <- row.names(df_all)
df_melt <- melt(df_all[,c("Gene", "avg_log2FC_moderate", "avg_log2FC_severe")])
df_melt$variable <- gsub(pattern = "avg_log2FC_", replacement = "", df_melt$variable)
df_melt$Gene <- factor(df_melt$Gene, levels = c(down_label_gene, up_label_gene))
df_melt$variable <- factor(df_melt$variable, levels = c("moderate", "severe"))

tmp <- df_all[, c("avg_log2FC_moderate", "avg_log2FC_severe")]
ord <- hclust( dist(tmp, method = "euclidean"), method = "ward.D" )$order
tmp$gene <- row.names(tmp)
tmp <- melt(tmp[rev(ord), ])
tmp$variable <- gsub(pattern = "avg_log2FC_", replacement = "", tmp$variable)
tmp$gene <- factor(tmp$gene, levels = rev(unique(tmp$gene)))
heatmap <- ggplot(tmp, aes(x = variable, y = gene, fill = value))+
  geom_tile(color = "black")+
    theme_classic() + 
      scale_fill_gradient2(midpoint = 0, low = "#2171b5", high = "#cb181d")+
      theme(text = element_text(size = 14, color = "black"),
            axis.title.y = element_blank(),
            axis.title.x = element_blank(),
            axis.text.x = element_text(size = 16, angle = 60, hjust = 1),
                plot.title = element_text(size = 16),
            strip.text.y =  element_text(angle = 0))+
      ggtitle("Avg logFC")

heatmap

4.4 Functional enrichment of core genes

plot_list <- list()

down <- df[df$avg_log2FC_moderate <(-0.1) & df$avg_log2FC_severe <(-0.1) & df$regulation %in% "DE_both", "Gene"]
up <- df[df$avg_log2FC_moderate >0.1 & df$avg_log2FC_severe >0.1 & df$regulation %in% "DE_both", "Gene"]

GSEA_core <- GSEA2(object = seurat_monocytes_late,
                 up_reg_genes = up, down_reg_genes = down,
                 name = "Dexa vs. Ctrl.",
                 GeneSets =c("GO", "Hallmark", "KEGG"),
                 GOntology = "BP",
                 pval.use = "p_val_adj",
                 pCorrection = "bonferroni", # choose the p-value adjustment method
                 pvalueCutoff = 0.05,
                 qvalueCutoff = 0.05 # set the q-value cutoff (FDR corrected)
)

GSEA_core$GOup$database <- "GO"
GSEA_core$HALLMARKup$database <- "HALLMARK"

df_up <- rbind(GSEA_core$GOup,
               GSEA_core$HALLMARKup)

df_up <- df_up[order(df_up$Count, decreasing = T),]
df_up$Description <- factor(df_up$Description, levels = rev(unique(df_up$Description)))
df_up$Description <- gsub(df_up$Description, pattern = "HALLMARK_", replacement = "")
df_up <- df_up[order(df_up$Count, decreasing = T),]
df_up$Description <- factor(df_up$Description, levels = rev(unique(df_up$Description)))

GSEA_core$GOdown$database <- "GO"
GSEA_core$HALLMARKdown$database <- "HALLMARK"

GSEA_core$GOdown$comparison <- NULL

df_down <- rbind(GSEA_core$GOdown,
               GSEA_core$HALLMARKdown)

df_down <- df_down[order(df_down$Count, decreasing = T),]
df_down$Description <- factor(df_down$Description, levels = rev(unique(df_down$Description)))

df_down$Description <- gsub(df_down$Description, pattern = "HALLMARK_", replacement = "")

df_down <- df_down[order(df_down$Count, decreasing = T),]
df_down$Description <- factor(df_down$Description, levels = rev(unique(df_down$Description)))

df_up$comparison <- NULL
df_up$reg <- "up"
df_down$reg <- "down"

df_both <- rbind(df_up[df_up$Count >0, ], head(df_down[order(df_down$Count, decreasing = T), ], 5))

df_both$reg <- factor(df_both$reg, levels = c("up", "down"))
df_both$Description <- factor(df_both$Description, levels = rev(unique(df_both$Description)) )

p <- ggplot(df_both, aes(x = reg, y = Description, fill = reg)) +
      geom_point(aes(size = Count), pch = 21, color = "black")+
             ggtitle(label = "functional Dexa core", subtitle =  "(p<0.05)")+
                theme_classic()+
  scale_fill_manual(values = c("#cb181d", "#2171b5"))+
                                            theme_bw()+
                                            theme_linedraw()+
                                            theme(panel.grid.major = element_line("white"),
                                                  panel.grid.minor = element_line("white"),
                                                  axis.title.y = element_text( size = 12),
                                                  axis.text.y = element_text( size = 12),
                                                  axis.title.x = element_text( size = 12),
                                                  axis.text.x = element_text( size = 12),
                                                  legend.text =  element_text( size = 12),
                                                  plot.title = element_text(face = "bold", size = 14))+
                                                  xlab("")+
                                                  ylab("Description")
p

4.5 Wang GC signature monocytes

enrichment of GC in-vitro signatures from Wang et al.

seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat_monocytes_late, group_severity %in% c("moderate", "severe"))

GC_wang <- read.csv("./files/Wang_et_al_Nature_GC_treatment_monocytes.csv", sep = ";")
GC_wang <- na.omit(GC_wang)
GC_wang_up <- GC_wang[order(GC_wang$Mo.DMSO_vs_TA_log2FoldChange, decreasing = T),][1:100,]
GC_wang_up <- GC_wang_up$gene
GC_wang_down <- GC_wang[order(GC_wang$Mo.DMSO_vs_TA_log2FoldChange, decreasing = F),][1:100,]
GC_wang_down <- GC_wang_down$gene

seurat <-  AddModuleScore(object = seurat,
                          features = list(GC_wang_up),
                          name = 'GC_wang_up')

seurat <-  AddModuleScore(object = seurat,
                          features = list(GC_wang_down),
                          name = 'GC_wang_down')

df <- seurat@meta.data[,c("donorID", "sampleID", "group", "GC_wang_up1")]
df %>% group_by(group) %>% summarize(mean = mean(GC_wang_up1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "group")

df <- df[order(df$scale_mean, decreasing = T),]
df$group <- factor(df$group, levels = c("moderate_ctrl", "moderate_Dexa", "severe_ctrl", "severe_Dexa"))

p1 <- ggplot(df, aes(x = group, y = GC_wang_up1))+
  geom_violin(aes(fill = scale_mean))+
  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
  scale_fill_gradient2(low = "blue", mid = "white",  high = "red", midpoint = 0)+
  theme_classic()+
  theme(panel.grid.major = element_line("white"),
        panel.grid.minor = element_line("white"),
        axis.title.y = element_text( size = 12),
        axis.text.y = element_text( size = 12),
        axis.title.x = element_text( size = 12),
        axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
        legend.text =  element_text( size = 12),
        plot.title = element_text(size = 14, hjust = 0.5))+
  xlab("")+
  ylab("Module Score")+
  ggtitle("GC_wang_up")+
  NoLegend()


df <- seurat@meta.data[,c("donorID", "sampleID", "group", "GC_wang_down1")]
df %>% group_by(group) %>% summarize(mean = mean(GC_wang_down1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "group")

df <- df[order(df$scale_mean, decreasing = T),]
df$group <- factor(df$group, levels = c("moderate_ctrl", "moderate_Dexa", "severe_ctrl", "severe_Dexa"))

p2 <- ggplot(df, aes(x = group, y = GC_wang_down1))+
  geom_violin(aes(fill = scale_mean))+
  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
  scale_fill_gradient2(low = "blue", mid = "white",  high = "red", midpoint = 0)+
  theme_classic()+
  theme(panel.grid.major = element_line("white"),
        panel.grid.minor = element_line("white"),
        axis.title.y = element_text( size = 12),
        axis.text.y = element_text( size = 12),
        axis.title.x = element_text( size = 12),
        axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
        legend.text =  element_text( size = 12),
        plot.title = element_text(size = 14, hjust = 0.5))+
  xlab("")+
  ylab("Module Score")+
  ggtitle("GC_wang_down")+
  NoLegend()

p1+p2

4.6 Monocyte UMAP

monocyte_levels <- c("IFIhi", "IL1Bhi", "S100Ahi", "CXCLhi",  "Dexa_response", "classical",  "THBDhi", "PF4+", "HBhi", "non_classical")
seurat_monocytes$cluster_labels_res.0.7 <- factor(seurat_monocytes$cluster_labels_res.0.7, levels = monocyte_levels)

DimPlot(seurat_monocytes, reduction = "umap", cols = colors_monocytes, group.by = "cluster_labels_res.0.7", order = F)

4.7 Monocyte state marker plot

seurat_monocytes$cluster_labels_res.0.7 <- factor(seurat_monocytes$cluster_labels_res.0.7, levels = monocyte_levels)

Idents(seurat_monocytes) <- "cluster_labels_res.0.7"
cluster.markers.wilcox_mono <- FindAllMarkers(object = seurat_monocytes, 
                                         only.pos = TRUE, 
                                         min.pct = 0.2, 
                                         logfc.threshold = 0.25,
                                         test.use = "wilcox", verbose = F
                                         )
top <- cluster.markers.wilcox_mono %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)

genes <- c("IFI6", "IFI27", "IFI44L", "XAF1", "MX1", 
           "IL1B", "CCL3", "CCL3L3", "CCL4", "CCL4L2", "TNF", 
           "S100A12", "S100A8", "S100A9", "RETN", "LGALS1",
           "CXCL1", "CXCL2", "CXCL3", "CXCL8", "MAFB", "HBEGF",
           "AREG", "IL1R2", "TSC22D3", "THBS1", "SAP30", "JDP2", "ZFP36L2", "EREG", "CXCR4", "CD163", "FKBP5", "KLF9", "PER1", 
           "HLA-DPA1", "HLA-DOB1", "HLA-DRA", "HLA-DRB1", "HLA-DQB1",
           "RGCC", "LMNA", "PHLDA1", "THBD", "CCL7", "EMP1", 
           "PPBP", "PF4", "TUBB1", "NRGN", "SPARC", 
           "HBB", "HBA1", "HBA2",
           "FCGR3A", "CX3CR1", "CDKN1C", "C1QA", "C1QB", "C1QC")

seurat_monocytes$cluster_labels_res.0.7 <- factor(seurat_monocytes$cluster_labels_res.0.7, levels = rev(monocyte_levels))
p <- DotPlot(seurat_monocytes,
        group.by = "cluster_labels_res.0.7", 
        features = genes)+ 
  RotatedAxis()+ 
  scale_color_gradientn(colors = rev(RColorBrewer::brewer.pal(n =20, name = "RdBu")))+
  theme(axis.title = element_blank()) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

p

4.8 Comparison to Schulte-Schrepping marker list

p <- DotPlot(seurat_monocytes,
        group.by = "cluster_labels_res.0.7", 
        features = c("MTND2P28", "HLA-DQB1", "HLA-DRB1", "HLA-DMB", "HLA-DPB1", "HLA-DRB5", # Classical Monocytes
                     "S100A12", "S100A8", "S100A9", "CXCL8", "ALOX5AP", "IRS2", "PLBD1",  # HLA-DRlo S100Ahi
                     "NRGN", "PF4", "PBPP", "MYL9", "CAVIN2", "SPARC", "TUBB1", "ITGA2B", # HLA-DRlo PF4+
                     "IFI6", "IFI27", "MX1",  "ISG15", "SELL", "CD163", "OAS2", "RETN",  # HLA-DRlo CD163hi
                     "HBB", "HBA2", "HBA1", "HBD",  # HLA-DRhi HBBhi
                     "RGCC",  "PHLDA1", "EMP1", "THBD", "SKG1", "ATP2B1", "CD83",   # HLA-DRhi CD83hi
                     "FCGR3A", "CDKN1C", "C1QB", "CX3CR1", "MS4A7", "C1QA", "C1QC" ))+ # non-classical Monocytes
  RotatedAxis()+ 
  scale_color_gradientn(colors = rev(RColorBrewer::brewer.pal(n =20, name = "RdBu")))+
  theme(axis.title = element_blank()) + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, face = "italic"))

p

4.9 Signature enrichment by monocyte state

response_glucocorticoid <- read.delim("./files/GO_term_summary_20230303_050304.txt", row.names = NULL)
response_glucocorticoid <- response_glucocorticoid$MGI.Gene.Marker.ID
response_glucocorticoid <- toupper(response_glucocorticoid)

response_steroid_hormone <- clusterProfiler::read.gmt("./files/GO_response_to_steroid_hormone.gmt")
response_steroid_hormone <- response_steroid_hormone$gene

seurat <- seurat_monocytes

seurat <-  AddModuleScore(object = seurat,
                          features = list(response_glucocorticoid),
                          name = 'response_glucocorticoid')
seurat <-  AddModuleScore(object = seurat,
                          features = list(response_steroid_hormone),
                          name = 'response_steroid')
seurat <-  AddModuleScore(object = seurat,
                          features = list(GC_wang_up),
                          name = 'GC_wang_up')
seurat <-  AddModuleScore(object = seurat,
                          features = list(GC_wang_down),
                          name = 'GC_wang_down')

df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "GC_wang_up1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(GC_wang_up1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")

df <- df[order(df$scale_mean, decreasing = T),]
levels_plot <- unique(df$cluster_labels_res.0.7)
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)

p1 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = GC_wang_up1))+
  geom_violin(aes(fill = scale_mean))+
  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
  scale_fill_gradient2(low = "blue", mid = "white",  high = "red", midpoint = 0)+
  theme_classic()+
  theme(panel.grid.major = element_line("white"),
        panel.grid.minor = element_line("white"),
        axis.title.y = element_text( size = 12),
        axis.text.y = element_text( size = 12),
        axis.title.x = element_text( size = 12),
        axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
        legend.text =  element_text( size = 12),
        plot.title = element_text(size = 14, hjust = 0.5))+
  xlab("")+
  ylab("Module Score")+
  ggtitle("GC_wang_up")+
  NoLegend()


df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "GC_wang_down1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(GC_wang_down1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")

df <- df[order(df$scale_mean, decreasing = T),]
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)

p2 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = GC_wang_down1))+
  geom_violin(aes(fill = scale_mean))+
  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
  scale_fill_gradient2(low = "blue", mid = "white",  high = "red", midpoint = 0)+
  theme_classic()+
  theme(panel.grid.major = element_line("white"),
        panel.grid.minor = element_line("white"),
        axis.title.y = element_text( size = 12),
        axis.text.y = element_text( size = 12),
        axis.title.x = element_text( size = 12),
        axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
        legend.text =  element_text( size = 12),
        plot.title = element_text(size = 14, hjust = 0.5))+
  xlab("")+
  ylab("Module Score")+
  ggtitle("GC_wang_down")+
  NoLegend()

df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "response_glucocorticoid1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(response_glucocorticoid1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)

p3 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = response_glucocorticoid1))+
  geom_violin(aes(fill = scale_mean))+
  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
  scale_fill_gradient2(low = "blue", mid = "white",  high = "red", midpoint = 0)+
  theme_classic()+
  theme(panel.grid.major = element_line("white"),
        panel.grid.minor = element_line("white"),
        axis.title.y = element_text( size = 12),
        axis.text.y = element_text( size = 12),
        axis.title.x = element_text( size = 12),
        axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
        legend.text =  element_text( size = 12),
        plot.title = element_text(size = 14, hjust = 0.5))+
  xlab("")+
  ylab("Module Score")+
  ggtitle("response_glucocorticoid")+
  NoLegend()


df <- seurat@meta.data[,c("donorID", "sampleID", "cluster_labels_res.0.7", "response_steroid1")]
df %>% group_by(cluster_labels_res.0.7) %>% summarize(mean = mean(response_steroid1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "cluster_labels_res.0.7")
df$cluster_labels_res.0.7 <- factor(df$cluster_labels_res.0.7, levels = levels_plot)

p4 <- ggplot(df, aes(x = cluster_labels_res.0.7, y = response_steroid1))+
  geom_violin(aes(fill = scale_mean))+
  stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
  scale_fill_gradient2(low = "blue", mid = "white",  high = "red", midpoint = 0)+
  theme_classic()+
  theme(panel.grid.major = element_line("white"),
        panel.grid.minor = element_line("white"),
        axis.title.y = element_text( size = 12),
        axis.text.y = element_text( size = 12),
        axis.title.x = element_text( size = 12),
        axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
        legend.text =  element_text( size = 12),
        plot.title = element_text(size = 14, hjust = 0.5))+
  xlab("")+
  ylab("Module Score")+
  ggtitle("response_steroid")+
  NoLegend()


p1+p2+p3+p4+plot_layout(ncol = 2)

5. Figure 3

5.1 Dexa response state distribution

seurat_monocytes_late <- subset(seurat_monocytes, order == "late")

seurat_quant <- subset(seurat_monocytes_late, group_severity %in% c("moderate", "severe"))

# quantify cells of each sample per cluster 
cells_cluster <- as.matrix(confusionMatrix(paste0(seurat_quant$cluster_labels_res.0.7),
                                           paste0(seurat_quant$sampleID)))
percentages <- round(t(t(cells_cluster)/colSums(cells_cluster))*100,3)

# transpose data for inspection
meta <- data.frame(seurat_quant@meta.data)
COVID_Dexa.percentages <- melt(percentages)
colnames(COVID_Dexa.percentages) <- c("cluster_labels_res.0.7","sampleID","percent")
idx <- match(COVID_Dexa.percentages$sampleID,meta$sampleID)
COVID_Dexa.percentages$group <- meta$group[idx]
COVID_Dexa.percentages$outcome <- meta$outcome[idx]

data <- COVID_Dexa.percentages[COVID_Dexa.percentages$cluster_labels_res.0.7 == "Dexa_response",]
data$group <- factor(data$group, levels = c("moderate_ctrl", "moderate_Dexa", 
                                             "severe_ctrl", "severe_Dexa")) 

p  <- ggplot(data, aes(fill = group, y=percent, x=group)) +
            geom_boxplot(outlier.size = 0, outlier.color = "white")+ 
            geom_point(pch = 21, color = "black", fill = "gray30" ,size=2,position = position_jitter(width = 0.1, height = 0))+
            theme_classic()+ 
            ggtitle("Dexa response state")+
            scale_fill_manual(values= c("#984EA3", "#5AAE61", "#984EA3", "#5AAE61")) +
            xlab("")+
            theme(axis.text.x = element_text(angle = 45, hjust =1, size = 12),
                  axis.text.y = element_text(size = 12), 
                  plot.title = element_text(size=16, face = "bold"))+
            ggpubr::stat_compare_means(comparisons = list(c("severe_ctrl", "severe_Dexa"), 
                                                          c("moderate_ctrl", "moderate_Dexa")), method = "wilcox.test")

p

5.2 Dexa response state by outcome

seurat_monocytes_late <- subset(seurat_monocytes, order == "late")

seurat_quant <- subset(seurat_monocytes_late, group %in% c("severe_Dexa"))

# quantify cells of each sample per cluster 
cells_cluster <- as.matrix(confusionMatrix(paste0(seurat_quant$cluster_labels_res.0.7),
                                           paste0(seurat_quant$sampleID)))
percentages <- round(t(t(cells_cluster)/colSums(cells_cluster))*100,3)

# transpose data for inspection
meta <- data.frame(seurat_quant@meta.data)
COVID_Dexa.percentages <- melt(percentages)
colnames(COVID_Dexa.percentages) <- c("cluster_labels_res.0.7","sampleID","percent")
idx <- match(COVID_Dexa.percentages$sampleID,meta$sampleID)
COVID_Dexa.percentages$group <- meta$group[idx]
COVID_Dexa.percentages$outcome <- meta$outcome[idx]

data <- COVID_Dexa.percentages[COVID_Dexa.percentages$cluster_labels_res.0.7 == "Dexa_response",]
data$outcome <- factor(data$outcome, levels = c("deceased", "survived")) 

p  <- ggplot(data, aes(fill = outcome, y=percent, x=outcome)) +
            geom_boxplot(outlier.size = 0, outlier.color = "white")+ 
            geom_point(pch = 21, color = "black", fill = "gray30" ,size=2,position = position_jitter(width = 0.1, height = 0))+
            theme_classic()+ 
            ggtitle("Dexa response state")+
            scale_fill_manual(values= c("#5B7674", "#57B391")) +
            xlab("")+
            theme(axis.text.x = element_text(angle = 45, hjust =1, size = 12),
                  axis.text.y = element_text(size = 12), 
                  plot.title = element_text(size=16, face = "bold"))+
            ggpubr::stat_compare_means(comparisons = list(c("deceased", "survived")), method = "wilcox.test")

p

## 5.3 Volcano severe survivors Dexa vs control

DEG_list <- list()

for(i in c("severe")){
  
  tmp <- subset(seurat_monocytes, subset = group_severity == i)
  tmp <- subset(tmp, subset = outcome == "survived")

  for (x in c("late")) {
    
    tmp_sample_order <- subset(tmp, subset = order == x)
    tmp_sample_order <- SetIdent(tmp_sample_order, value = "group")
    
      ctrl <- unique(tmp_sample_order$group)[grep(unique(tmp_sample_order$group), pattern = "ctrl")]
      
      idx <- ifelse(grep(unique(tmp_sample_order$group), pattern = "ctrl") == 1, 2, 1)
      treat <- unique(tmp_sample_order$group)[idx]
      
      tryCatch({
        DEG<-FindMarkers(tmp_sample_order,
                         ident.1 =  treat,
                         ident.2 = ctrl,
                         min.pct = 0.1,
                         logfc.threshold = 0,
                         test.use = "wilcox",
                         slot = "data",
                         only.pos = F,
                         verbose = F)
        DEG$Gene <- rownames(DEG)
        DEG$group_severity <- paste(i)
        DEG$sample_comparison <- paste(x)
        DEG$comparison <- paste0(i, "_", x, "_treatment_vs_ctrl")
        DEG_list[[paste0(i, "_", x, "_treatment_vs_ctrl")]]<-DEG
      }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
 }
}

up_reg <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
                                                         DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC>0,]
up_reg <- up_reg[order(up_reg$avg_log2FC, -up_reg$p_val_adj, decreasing = T),]$Gene
            
down_reg <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
                                                           DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$avg_log2FC, -down_reg$p_val_adj, decreasing = F),]$Gene
              
              
up_reg_pval <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
                                                              DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC>0,]
up_reg_pval <- up_reg_pval[order(up_reg_pval$p_val_adj, decreasing = F),]$Gene
              
down_reg <- DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
                                                           DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$p_val_adj, decreasing = F),]$Gene
              

p <-ggplot(DEG_list$severe_late_treatment_vs_ctrl,aes(x=avg_log2FC,y=-log10(p_val_adj))) + 
          geom_point(pch = 21, fill = "darkgrey", color = "black", size = 2) +
          geom_text_repel(data=DEG_list$severe_late_treatment_vs_ctrl[c("CST3", "CD74", "HLA-DRB1", "HLA-DRA", 
                                                                        "AREG", "HLA-DPA1", "TSC22D3",
                                                                        "THBS1", "CEBPD",  "HLA-DPB1", 
                                                                        "IL1R2", "MT-ND3", "HLA-F", "CXCR4", 
                                                                        "ZFP36L2"),],
                          max.overlaps = 20, 
                          aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),
                          size=3, color = "black")+
         geom_text_repel(data=DEG_list$severe_late_treatment_vs_ctrl[c("CCL4","S100A9", "S100A8", "S100A12", 
                                                                       "CCL3", "CCL4L2", "CCL3L3", "IL1B", 
                                                                       "NAMPTP1","IL1RN", "RETN", "CLU",
                                                                       "CXCL3", "G0S2", "CXCL2", "PHLDA1",
                                                                       "RGCC", "CXCL8", "HP", "BCL2A1", "CD36",
                                                                       "FTH1", "CST7"),],
                         max.overlaps = 20, 
                         aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),size=3, 
                         color = "black")+
        geom_point(data=DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
                                                                               DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC>0.2,],
                                aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#57B391", size = 2, pch = 21, color = "black")+
        geom_point(data=DEG_list$severe_late_treatment_vs_ctrl[DEG_list$severe_late_treatment_vs_ctrl$p_val_adj<.05&
                                                                               DEG_list$severe_late_treatment_vs_ctrl$avg_log2FC<c(-0.2),],
                                aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#F08FB6", size = 2, pch = 21, color = "black")+
        ggtitle(paste0("survivor DEG analysis (severe COVID-19)"))+
        theme_classic()+
        theme_bw()+ 
        theme_linedraw()+
        theme(panel.grid.major = element_line("white"), 
              panel.grid.minor = element_line("white"),
              axis.title.y = element_text( size = 8),
              axis.text.y = element_text( size = 8),
              axis.title.x = element_text( size = 8),
              axis.text.x = element_text( size = 8),
              legend.title = element_blank(),
              legend.position = "none", 
              plot.title = element_text(face = "bold", size = 10))+
        xlab("log2FC")+
        ylab("-log10(p_val_adj)")+
        geom_hline(yintercept = 0)+ 
        geom_vline(xintercept = 0)+
        xlim(c(-2, 2.5))+
        ylim(c(0, 220))
            
p   

5.4 volcano of selected DEG Dexa survivors

seurat <- subset(seurat_monocytes, order == "late" & outcome == "survived" & group_severity %in% c("severe"))

p1 <- VlnPlot(seurat, features = "HLA-DRB1", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p2 <- VlnPlot(seurat, features = "HLA-DRA", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p3 <- VlnPlot(seurat, features = "HLA-DPA1", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p4 <- VlnPlot(seurat, features = "CD74", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p5 <- VlnPlot(seurat, features = "S100A8", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p6 <- VlnPlot(seurat, features = "S100A9", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p7 <- VlnPlot(seurat, features = "S100A12", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p8 <- VlnPlot(seurat, features = "CCL4", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p9 <- VlnPlot(seurat, features = "IL1B", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p10 <- VlnPlot(seurat, features = "CXCL3", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p11 <- VlnPlot(seurat, features = "CXCL2", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())
p12 <- VlnPlot(seurat, features = "CXCL8", group.by = "group", pt.size = 0, cols = c("#F08FB6", "#57B391"))+ 
              stat_summary(fun.y = median, geom='point', size = 10, colour = "black", shape = 95)+
              NoLegend()+
              theme(axis.title.y = element_blank(), axis.title.x = element_blank())


p1+p2 + p5+p6 + p9+p10 +  
p3+p4 + p7+p8 + p11+p12 + plot_layout(nrow = 2, ncol = 6)

5.5 S100A state distribution

seurat_monocytes_late <- subset(seurat_monocytes, order == "late")

seurat_quant <- subset(seurat_monocytes_late, group %in% c("severe_Dexa"))

# quantify cells of each sample per cluster 
cells_cluster <- as.matrix(confusionMatrix(paste0(seurat_quant$cluster_labels_res.0.7),
                                           paste0(seurat_quant$sampleID)))
percentages <- round(t(t(cells_cluster)/colSums(cells_cluster))*100,3)

# transpose data for inspection
meta <- data.frame(seurat_quant@meta.data)
COVID_Dexa.percentages <- melt(percentages)
colnames(COVID_Dexa.percentages) <- c("cluster_labels_res.0.7","sampleID","percent")
idx <- match(COVID_Dexa.percentages$sampleID,meta$sampleID)
COVID_Dexa.percentages$group <- meta$group[idx]
COVID_Dexa.percentages$outcome <- meta$outcome[idx]

data <- COVID_Dexa.percentages[COVID_Dexa.percentages$cluster_labels_res.0.7 == "S100Ahi",]
data$outcome <- factor(data$outcome, levels = c("deceased", "survived")) 

p  <- ggplot(data, aes(fill = outcome, y=percent, x=outcome)) +
            geom_boxplot(outlier.size = 0, outlier.color = "white")+ 
            geom_point(pch = 21, color = "black", fill = "gray30" ,size=2,position = position_jitter(width = 0.1, height = 0))+
            theme_classic()+ 
            ggtitle("S100Ahi state")+
            scale_fill_manual(values= c("#5B7674", "#57B391")) +
            xlab("")+
            theme(axis.text.x = element_text(angle = 45, hjust =1, size = 12),
                  axis.text.y = element_text(size = 12), 
                  plot.title = element_text(size=16, face = "bold"))+
            ggpubr::stat_compare_means(comparisons = list(c("deceased", "survived")), method = "wilcox.test")

p

5.6 Dexa outcome DEG volcano (severe COVID-19)

seurat <- subset(seurat_monocytes, subset = group_severity == "severe" & order == "late")
seurat <- SetIdent(seurat, value = "group_outcome")

DEG_deceased <- list()

x= "late"
tryCatch({
  DEG<-FindMarkers(seurat,
                   ident.1 =  "severe_Dexa_survived",
                   ident.2 = "severe_Dexa_deceased",
                   min.pct = 0.1,
                   logfc.threshold = 0,
                   test.use = "wilcox",
                   slot = "data",
                   only.pos = F,
                   verbose = F)
  DEG$Gene <- rownames(DEG)
  DEG_deceased[[ paste0(x, "_treatment_alive_vs_deceased")]]<-DEG
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})

df <- DEG_deceased$late_treatment_alive_vs_deceased
    
up_reg <- df[df$p_val_adj<.05&df$avg_log2FC>0,]
up_reg <- up_reg[order(up_reg$avg_log2FC, -up_reg$p_val_adj, decreasing = T),]$Gene
down_reg <- df[df$p_val_adj<.05&df$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$avg_log2FC, -down_reg$p_val_adj, decreasing = F),]$Gene
up_reg_pval <- df[df$p_val_adj<.05&df$avg_log2FC>0,]
up_reg_pval <- up_reg_pval[order(up_reg_pval$p_val_adj, decreasing = F),]$Gene
down_reg <- df[df$p_val_adj<.05&df$avg_log2FC<0,]
down_reg <- down_reg[order(down_reg$p_val_adj, decreasing = F),]$Gene
              
p <- ggplot(df,aes(x=avg_log2FC,y=-log10(p_val_adj))) + 
            geom_point(pch = 21, fill = "darkgrey", color = "black", size = 2) +
            geom_text_repel(data=df[c("S100A8", "S100A9", "S100A12", "PKM", "CLU", 
                                      "LGALS1", "MCEMP1", "RETN", "MT-ATP8", 
                                      "TMEM176B", "PLAC8", "SERPINB2"),],
                                max.overlaps = 20, 
                                aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),size=3, color = "black")+
              
             geom_text_repel(data=df[c("HLA-DRB1", "CD74", "HLA-DRA", 
                                       "HLA-DPA1", "CEBPD", "HLA-DPB1", 
                                       "ZFP36L2", "HLA-DQB1", "IFI6", 
                                       "AREG", "CCL3","CCL4", "CCL3L3",  
                                       "IL1R2",    "HLA-DMB",  "HLA-DRB5", 
                                       "RPL34", "CST3", "RPS13"),],
                                max.overlaps = 20,
                                aes(x=avg_log2FC,y=-log10(p_val_adj),label=Gene),size=3, color = "black")+
             geom_point(data=df[df$p_val_adj<.05&df$avg_log2FC>0.2,],
                                aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#57B391", size = 2, pch = 21, color = "black")+
             geom_point(data=df[df$p_val_adj<.05&df$avg_log2FC<c(-0.2),],
                                aes(x=avg_log2FC,y=-log10(p_val_adj)), fill = "#5B7674", size = 2, pch = 21, color = "black")+
             ggtitle(paste0("Dexa outcome DEG analysis (severe COVID-19)"))+
             theme_classic()+
             theme_bw()+ 
             theme_linedraw()+
             theme(panel.grid.major = element_line("white"),
                  panel.grid.minor = element_line("white"),
                  axis.title.y = element_text( size = 8),
                  axis.text.y = element_text( size = 8),
                  axis.title.x = element_text( size = 8),
                  axis.text.x = element_text( size = 8),
                  legend.title = element_blank(),
                  legend.position = "none", 
                  plot.title = element_text(face = "bold", size = 10))+
                  xlab("log2FC")+
                  ylab("-log10(p_val_adj)")+
             geom_hline(yintercept = 0)+ 
             geom_vline(xintercept = 0)+
             xlim(c(-1.3, 1.6))+
             ylim(c(0, 280))
 
p

5.7 selected outcome DEG

seurat_monocytes_late <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat_monocytes_late, group_severity %in% "severe")
seurat <- subset(seurat, group_outcome %in% "severe_ctrl_deceased", invert = T)

seurat$group_outcome <- factor(seurat$group_outcome, levels = rev(c("severe_ctrl_survived", 
                                                                    "severe_Dexa_deceased", 
                                                                    "severe_Dexa_survived")))

p <- DotPlot(seurat, group.by = "group_outcome", 
             features = c("S100A8", "S100A9", "S100A12",
                          "CD74", "HLA-DRB1", "HLA-DRA", 
                          "HLA-DPA1", "HLA-DPB1", "HLA-DQB1",
                          "HLA-DMB", "ZFP36L2", "IFI6", "AREG"), 
             col.max = 1.5, col.min = -1.5)+ 
             theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
             scale_color_gradient2(low = "blue", mid = "white", high = "red")

p

check for correlation of the selected genes.

seurat$group_outcome_donor <- paste(seurat$group_outcome, seurat$donorID, sep = "_")

df <- AverageExpression(seurat, features = c("S100A8", "S100A9", "S100A12", 
                                                                        "CD74", "HLA-DRB1", "HLA-DRA", "HLA-DPA1", "HLA-DPB1", "HLA-DQB1", "HLA-DMB", "ZFP36L2", "IFI6", "AREG"
                                                                        ), group.by = "group_outcome")$RNA
df
##          severe_Dexa_survived severe_Dexa_deceased severe_ctrl_survived
## S100A8             119.475170          206.7166969          203.6020996
## S100A9             103.012413          158.3877037          172.6169736
## S100A12             11.206291           22.2545702           17.8782846
## CD74                45.761604           21.9769300           27.6624037
## HLA-DRB1            12.116669            3.6801748            5.7532979
## HLA-DRA             17.651813            7.1160404           10.2649664
## HLA-DPA1             4.711761            1.3969531            2.4362415
## HLA-DPB1             2.601729            0.7854829            1.2535002
## HLA-DQB1             1.841001            0.5625530            1.5580807
## HLA-DMB              1.578108            0.6232262            0.8807122
## ZFP36L2              4.901420            2.1062185            2.2731095
## IFI6                 6.430200            2.1234939            3.8872294
## AREG                 4.219793            1.2213803            2.0480760
df <- as.data.frame(t(scale(t(df))))

df <- as.data.frame(df)

cor_df <- cor(df)
cor_df <- cor_df[c("severe_ctrl_survived", "severe_Dexa_deceased", "severe_Dexa_survived"), ]

pheatmap(cor_df[,"severe_ctrl_survived"], 
         cluster_cols = F, 
         cluster_rows = F, 
         display_numbers = T, 
         number_color = "black")

5.8 HLA-DRlo S100Ahi Schulte-Schrepping signature enrichment

seurat_Schulte <- readRDS("./files/seurat_COVID19_monocytes_jonas_FG_2020-07-29.rds")
seurat_Schulte@meta.data$cell_ID <- row.names(seurat_Schulte@meta.data)
seurat_Schulte$dataset <- "Schulte"
seurat_Schulte$cluster_labels_res.0.2<- seurat_Schulte$RNA_snn_res.0.2
Idents(seurat_Schulte) <- "cluster_labels_res.0.2"
seurat_Schulte <- RenameIdents(object = seurat_Schulte, 
                             `0` = "HLA-DRlo_S100Ahi", 
                             `1` = "HLA-DRhi_CD83hi", 
                             `2` = "HLA-DRlo_CD163hi",
                             `3` = "Classical_Monocytes",
                             `4` = "Non-classical_Monocytes",
                             `5` = "HLA-DRhi_HBBhi",
                             `6` = "HLA-DRlo_PF4+")
seurat_Schulte@meta.data$cluster_labels_res.0.2 <- Idents(seurat_Schulte)

marker <- FindMarkers(seurat_Schulte, ident.1 = "HLA-DRlo_S100Ahi", 
                      group.by = "cluster_labels_res.0.2", 
                      logfc.threshold = 0.25, min.pct = 0.2, only.pos = T)
marker <- marker[marker$p_val_adj < 0.05, ]
marker <- marker[order(marker$avg_log2FC, decreasing = T), ]
Schulte_S100_top50 <- head(row.names(marker), 50)

seurat <- subset(seurat_monocytes, order == "late")
seurat <- subset(seurat, group_severity %in% c("severe"))

seurat <-  AddModuleScore(object = seurat,
                          features = list(Schulte_S100_top50),
                          name = 'HLA_DRlo_S100Ahi_signature')

df <- seurat@meta.data[,c("donorID", "sampleID", "group_outcome", "HLA_DRlo_S100Ahi_signature1")]
df %>% group_by(group_outcome) %>% summarize(mean = mean(HLA_DRlo_S100Ahi_signature1)) -> df.mean
df.mean$scale_mean <- scale(df.mean$mean)
df <- merge(df.mean, df, by = "group_outcome")
df <- df[order(df$scale_mean, decreasing = T),]
df$group_outcome <- factor(df$group_outcome,  levels = c("severe_ctrl_deceased",
                                                         "severe_ctrl_survived",
                                                         "severe_Dexa_deceased",
                                                         "severe_Dexa_survived"))

p <- ggplot(df, aes(x = group_outcome, y = HLA_DRlo_S100Ahi_signature1))+
            geom_violin(aes(fill = scale_mean))+
            stat_summary(fun.y = mean, geom='point', size = 10, colour = "black", shape = 95)+
            scale_fill_gradient2(low = "blue", mid = "white",  high = "red", midpoint = -0.5)+
            theme_classic()+
            theme(panel.grid.major = element_line("white"),
                  panel.grid.minor = element_line("white"),
                  axis.title.y = element_text( size = 12),
                  axis.text.y = element_text( size = 12),
                  axis.title.x = element_text( size = 12),
                  axis.text.x = element_text( size = 12, angle = 45, hjust = 1),
                  legend.text =  element_text( size = 12),
                  plot.title = element_text(size = 10, hjust = 0.5))+
            xlab("")+
            ylab("Module Score")+
            ggtitle("HLA-DRlo S100Ahi signature (Schulte-Schrepping, n=50)")

p

6. Enviornment info

R.version
##                _                           
## platform       x86_64-pc-linux-gnu         
## arch           x86_64                      
## os             linux-gnu                   
## system         x86_64, linux-gnu           
## status                                     
## major          4                           
## minor          1.0                         
## year           2021                        
## month          05                          
## day            18                          
## svn rev        80317                       
## language       R                           
## version.string R version 4.1.0 (2021-05-18)
## nickname       Camp Pontanezen
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.8.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] fmsb_0.7.5            viridis_0.6.1         viridisLite_0.4.0    
##  [4] cowplot_1.1.1         gridExtra_2.3         pheatmap_1.0.12      
##  [7] patchwork_1.1.1       circlize_0.4.13       tidyr_1.1.4          
## [10] org.Hs.eg.db_3.13.0   AnnotationDbi_1.54.1  IRanges_2.26.0       
## [13] S4Vectors_0.30.2      Biobase_2.52.0        BiocGenerics_0.38.0  
## [16] RColorBrewer_1.1-2    useful_1.2.6          Matrix_1.5-4.1       
## [19] ggnewscale_0.4.8      ggrepel_0.9.1         ggpubr_0.4.0         
## [22] ComplexHeatmap_2.8.0  clusterProfiler_4.0.5 harmony_0.1.0        
## [25] Rcpp_1.0.7            dplyr_1.0.7           ggplot2_3.3.5        
## [28] tibble_3.1.5          SeuratObject_4.1.3    Seurat_4.3.0         
## [31] stringr_1.4.0         reshape2_1.4.4       
## 
## loaded via a namespace (and not attached):
##   [1] scattermore_0.7        bit64_4.0.5            knitr_1.33            
##   [4] irlba_2.3.3            data.table_1.14.2      KEGGREST_1.32.0       
##   [7] RCurl_1.98-1.5         doParallel_1.0.16      generics_0.1.0        
##  [10] RSQLite_2.2.8          shadowtext_0.0.9       RANN_2.6.1            
##  [13] future_1.12.0          bit_4.0.4              enrichplot_1.12.3     
##  [16] spatstat.data_3.0-1    httpuv_1.6.3           assertthat_0.2.1      
##  [19] xfun_0.24              hms_1.1.1              jquerylib_0.1.4       
##  [22] evaluate_0.14          promises_1.2.0.1       fansi_0.5.0           
##  [25] readxl_1.3.1           igraph_1.2.6           DBI_1.1.1             
##  [28] htmlwidgets_1.5.4      spatstat.geom_3.2-1    purrr_0.3.4           
##  [31] ellipsis_0.3.2         backports_1.2.1        deldir_1.0-5          
##  [34] vctrs_0.3.8            Cairo_1.5-12.2         ROCR_1.0-11           
##  [37] abind_1.4-5            cachem_1.0.6           withr_2.4.2           
##  [40] ggforce_0.3.3          progressr_0.13.0       sctransform_0.3.5     
##  [43] treeio_1.16.2          goftest_1.2-3          cluster_2.1.2         
##  [46] DOSE_3.18.3            ape_5.5                lazyeval_0.2.2        
##  [49] crayon_1.4.1           spatstat.explore_3.2-1 pkgconfig_2.0.3       
##  [52] labeling_0.4.2         tweenr_1.0.2           GenomeInfoDb_1.28.4   
##  [55] vipor_0.4.5            nlme_3.1-152           rlang_1.1.1           
##  [58] globals_0.14.0         lifecycle_1.0.1        miniUI_0.1.1.1        
##  [61] downloader_0.4         ggrastr_0.2.3          cellranger_1.1.0      
##  [64] polyclip_1.10-0        matrixStats_0.61.0     lmtest_0.9-38         
##  [67] aplot_0.1.1            carData_3.0-4          zoo_1.8-9             
##  [70] beeswarm_0.4.0         ggridges_0.5.3         GlobalOptions_0.1.2   
##  [73] png_0.1-7              rjson_0.2.20           bitops_1.0-7          
##  [76] KernSmooth_2.23-20     Biostrings_2.60.2      blob_1.2.2            
##  [79] shape_1.4.6            qvalue_2.24.0          spatstat.random_3.1-5 
##  [82] rstatix_0.7.0          gridGraphics_0.5-1     ggsignif_0.6.3        
##  [85] scales_1.1.1           memoise_2.0.0          magrittr_2.0.1        
##  [88] plyr_1.8.6             ica_1.0-2              zlibbioc_1.38.0       
##  [91] compiler_4.1.0         scatterpie_0.1.7       clue_0.3-60           
##  [94] fitdistrplus_1.1-6     cli_3.6.1              XVector_0.32.0        
##  [97] listenv_0.8.0          pbapply_1.5-0          MASS_7.3-54           
## [100] tidyselect_1.1.1       stringi_1.7.5          forcats_0.5.1         
## [103] highr_0.9              yaml_2.2.1             GOSemSim_2.18.1       
## [106] sass_0.4.0             fastmatch_1.1-3        tools_4.1.0           
## [109] future.apply_1.2.0     rio_0.5.27             rstudioapi_0.13       
## [112] foreach_1.5.1          foreign_0.8-81         farver_2.1.0          
## [115] Rtsne_0.15             ggraph_2.0.5           digest_0.6.28         
## [118] shiny_1.7.1            car_3.0-11             broom_0.7.9           
## [121] later_1.3.0            RcppAnnoy_0.0.19       httr_1.4.2            
## [124] colorspace_2.0-2       tensor_1.5             reticulate_1.22       
## [127] splines_4.1.0          uwot_0.1.14            yulab.utils_0.0.4     
## [130] tidytree_0.3.5         spatstat.utils_3.0-3   graphlayouts_0.7.1    
## [133] sp_1.6-0               ggplotify_0.1.0        plotly_4.10.0         
## [136] xtable_1.8-4           jsonlite_1.7.2         ggtree_3.0.4          
## [139] tidygraph_1.2.0        ggfun_0.0.4            R6_2.5.1              
## [142] pillar_1.6.3           htmltools_0.5.2        mime_0.12             
## [145] glue_1.4.2             fastmap_1.1.0          BiocParallel_1.26.2   
## [148] codetools_0.2-18       fgsea_1.18.0           utf8_1.2.2            
## [151] lattice_0.20-44        bslib_0.3.1            spatstat.sparse_3.0-1 
## [154] ggbeeswarm_0.6.0       curl_4.3.2             leiden_0.3.9          
## [157] zip_2.2.0              GO.db_3.13.0           openxlsx_4.2.4        
## [160] survival_3.2-11        limma_3.48.3           rmarkdown_2.11        
## [163] munsell_0.5.0          DO.db_2.9              GetoptLong_1.0.5      
## [166] GenomeInfoDbData_1.2.6 iterators_1.0.13       haven_2.4.3           
## [169] gtable_0.3.0